Systematic investigation of a family of gradient-dependent functionals for solids 
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Eleven density functionals are compared with regard to their performance for the lattice constants 
of solids. We consider standard functionals, such as the local-density approximation and the Perdew- 
Burke-Ernzerhof (PBE) generalized-gradient approximation (GGA), as well as variations of PBE 
GGA, such as PBEsol and similar functionals, PBE-type functionals employing a tighter Lieb- 
Oxford bound, and combinations thereof. Several of these variations are proposed here for the first 
time. On a test set of 60 solids we perform a system-by-system analysis for selected functionals 
and a full statistical analysis for all of them. The impact of restoring the gradient expansion and of 
tightening the Lieb-Oxford bound is discussed, and confronted with previous results obtained from 
other codes, functionals or test sets. No functional is uniformly good for all investigated systems, 
but surprisingly, and pleasingly, the simplest possible modifications to PBE turn out to have the 
most beneficial effect on its performance. The atomization energy of molecules was also considered 
and on a testing set of six molecules, we found that the PBE functional is clearly the best, the 
others leading to strong overbinding. 

PACS numbers: 71.15.Mb, 71.15.Nc, 31.15.E- 
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I. INTRODUCTION 

Modern electronic-structure theory^ relies to a very 
large extent on density- functional theory (DFT)^- The 
utility of DFT, in turn, depends crucially on the availabil- 
ity of approximations to the exchange-correlation (xc) 
functional that are sufficiently reliable and sufficiently 
simple to implement i - — 

As a consequence, a large number of approximate xc 
functionals have been developed. Only a few of these, 
however, have found widespread application, and essen- 
tially just two of them account for the large majority 
of applications of DFT in solid-state physics: the local- 
density approximation (LDA) and the Perdew-Burke- 
Ernzerhof (PBE) form of the generalized-gradient ap- 
proximation (GGA)ji 

Among the main problems of these functionals is that 
lattice constants are systematically and consistently un- 
derestimated by LDA and overestimated by PBE. LDA 
lattice constants are typically about 1 — 5% too short, 
while PBE lattice constants are too long by almost the 
same margin. Many other quantities, such as the unit-cell 
geometry and volume, the cohesive energy, bulk mod- 
ulus, compressibility, phonon frequencies, sound veloc- 



ity, elastic constants, Debye temperature, the pressure- 
dependence of all these quantities, surface reconstruction 
energies, the possibility of structural phase transitions, 
etc., depend crucially on the lattice constant. There- 
fore, the difficulty of LDA and PBE in predicting quan- 
titatively reliable lattice constants is a crucial problem 
standing in the way of further applications of DFT to 
solids. 

Until quite recently, no generally applicable solution 
to this problem was in sight, and the very voluminous 
literature on better xc functionals (e.g., the hybrid func- 
tionals) largely focused on finite systems (see Refs. [H 
and |9| for recent reviews). For solids, however, these 
functionals do not perform that well in every situation 
and/or lead to very expensive calculations. For instance, 
the very popular hybrid functional B3LY P 10 i n is rather 
hard to implement for solids, in particular for metals, 
and the effort does not seem to pay off, as resulting 
lattice constants overestimate experimental values by 
about as much as PBE. 12 Similarly, semi-empiricali^ and 
nonempiricalii meta-GGA functionals (slightly more ex- 
pensive than GGAs) produce littl o 14 ! 15 or no^ improve- 
ment for lattice constants. Although many other func- 
tionals have been tried over the years, LDA and PBE 
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remained, until very recently, the de facto standard DFT 
approach for the determination of structural properties 
of solids and nanostructures. 

Recently, however, the field of functional construction 
for solids has gained new impetus through the develop- 
ment of AM05,— ~— a radically new type of density func- 
tional based on the subsystem approach and the Airy 
gas, and the Wu-Cohen (WC) GGA ) 19 i 20 which employs 
a simple but efficient modification of the PBE exchange 
enhancement factor that makes it more reliable for solid- 
state properties. An even simpler modification of PBE 
is PBEsol^I which differs from original PBE only in the 
values of two parameters. 

These developments have rekindled the interest in the 
development of better density functionals for solids. Sev- 
eral such recently developed functionals, AM05, WC, 
PBEsol, and the second-order GGA (SOGGA) of Zhao 
and Truhlar— have been systematically tested and com- 
pared to LDA, PBE and TPSS meta-GGA in a previ- 
ous publication by three of us (PH, FT, and PB)^ A 
large test set (60 solids) and a very accurate all-electron 
implementation of the Kohn-Sham (KS) equations (the 
WIEN2k code2£) allowed a detailed investigation of the 
performance of each of these functionals. Overall, no 
clear winner has emerged from the comparison, but the 
new GGA functionals improve over LDA and PBE for 
many solids and give smaller mean errors. In Ref. [HI, 
we reported a detailed analysis of the functionals, which 
shed light on some of the trends observed in the lattice 
constants. 

In an independent work, three of us (LSP, AJRdS, and 
KC) noted that the step that led from PBE to PBEsol is 
not unique, and allows several variations. 2 ~ In fact, PBE 
and PBEsol turned out to be just two particular mem- 
bers of a family of functionals each of which takes its 
parameters, /3 and /i, from a different constraint. The re- 
sulting two-parameter family of functionals, collectively 
denoted PBE(/3, p,), has been tested for atoms, molecules 
and solids in Ref. [25|. The calculations for solids per- 
formed in that work employed pseudopotentials, which is 
the standard approach for very large systems with many 
inequivalent sites, but introduces an additional source 
of errors and complicates an unbiased assessment of the 
performance of each functional. 

In still other work, two of us (MMO and KC) initiated 
an investigation of the Lieb-Oxford (LO) bound ) 26 i 27 
a fundamental property of the quantum mechanics of 
Coulomb-interacting systems according to which the ex- 
act xc energy is bounded from below by a simple lo- 
cal density functional that is proportional to the LDA 
for exchange. An estimate of the proportionality factor, 
A, is a parameter in several modern density function- 
als, among them SOGGA, TPSS, WC, as well as PBE, 
PBEsol and all other members of the PBE(/3, /x) fam- 
ily. Numerical and analytical investigations^ 7 - - — strongly 
suggest that the value adopted in standard density func- 
tionals Alo = 2.273 is too large and should be replaced 
by Ael = 1-9555. Of the functionals listed above only 



SOGGA makes use of this tighter bound. Consequences 
of a tighter LO bound in PBE calculations for molecular 
systems have been explored in Ref. Hj| but that work did 
not consider solids and did not include variations in /3 
and fi. We also mention the work of Peltzcr y Blanca 
et al£^ who concluded that reducing the value of A in 
PBE leads to better results for the equilibrium volume of 
Ad and 5d transition metals as also shown in the present 
work. 

In the present paper we now bring all these devel- 
opments together. We propose and study the three- 
parameter family of density functionals PBE(/3, /i, A), ex- 
plore all meaningful nonempirical combinations of these 
parameters that are available, implement the resulting 
ten functionals in the WIEN2k all-electron code,— and 
test them on the large set of 60 solids from Ref. O In 
addition we implemented these functionals in the deMon 
code22 and tested atomization energies on a small but 
representative set of six molecules.— 

This paper is organized as follows. In Sec. [TT] we de- 
scribe the PBE(/3,/i, A) family of functionals, indicating 
the possible values and sources of each of its three pa- 
rameters. Section [TTT] is devoted to a system-by-system 
comparison of three members of the family that differ in 
just one constraint from original PBE. In Sec. lIVI we then 
present results from a statistical analysis of the full set 
of ten PBE-type functionals and LDA, for all 60 solids. 
This section also contains a comparison of our results 
with those from several other published tests of similar 
functionals, among them various using different codes 
and different test sets. Section [V] analyses our results 
separately for elements and compounds (metallic tran- 
sition metal compounds, semiconductors and insulators) 
and Sec. IVII reports the performance of the PBE(/3, /i, A) 
functionals for atomization energies of small molecules. 
Finally Sec. IVIII contains our conclusions. 



II. THE PBE(,fl,^A) FAMILY OF 
FUNCTIONALS 

The structure of PBE is explained in the original 
reference^ and more details are given in the review 
literature^ In the interest of conciseness, we thus refrain 
from repeating the explicit expression of this widely used 
functional and directly focus on its parameters and their 
possible modifications. 

PBE contains nonempirical parameters, whose numer- 
ical values are obtained by requiring that the functional 
obeys known universal constraints. Two of them, k and 
/j,, appear in the exchange functional, E x , and one, /3, 
appears in the correlation functional, E c . 

In the original construction of PBE,— the parameter f3 
is chosen such that in the high-density limit E^ BE re- 
covers the second-order gradient expansion of the corre- 
lation energy of spatially weakly varying systems. The 
requirement that the combined xc functional reproduces 
the LDA jellium response function (which is accurate) 
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implies 

M=y/?, (1) 

which fixes /i. The third parameter, k, was deter- 
mined such that E X BE alone obeys the Lieb-Oxford lower 
bound 26 , on the xc energy. This implies 

« = 2^3 " 1 = °- 804 > ( 2 ) 

where Alo = 2.273 is an estimate of the Lieb-Oxford 
constant A obtained in Ref. l26l 

This particular choice of constraints proved to be 
enormously successful, and PBE is one of the most 
widely used density functionals across physics and chem- 
istry. Nevertheless, the choice is clearly not unique, 
and has recently been reconsidered along two indepen- 
dent lines. To discuss these, we introduce the notation 
PBE(/3, jit, A), where the parameters can be replaced ei- 
ther by their numerical values or by symbols indicating 
the source of these values. Hence, original PBE becomes 
PBE(G C , J r , LO) indicating that (3 comes from the gradi- 
ent expansion of E c , /x from the jellium response function, 
while for A the original Lieb-Oxford estimate is adopted. 

The first line of thought originates with the PBEsol 
functional, designed specifically to improve on PBE for 
solids. 21 To construct PBEsol it was argued that for 
solids the gradient expansion of the exchange functional 
is expected to be more important than that of the corre- 
lation functional. Consequently /i, which appears in the 
exchange energy, is chosen in PBEsol such as to repro- 
duce the second-order gradient expansion of E x . The pa- 
rameter (3, appearing in the correlation energy, is deter- 
mined in PBEsol by requiring that jellium surface ener- 
gies are accurately reproduced. In our notation, PBEsol 
becomes PBE( J a , G x , LO). PBEsol has been extensively 
teste d 15 ' 21 ' 35 " — and the results have vindicated the re- 
vised choice of constraints, as PBEsol indeed provides 
significant improvement on PBE for solids (at the ex- 
pense of worsening the results for smaller molecular sys- 
tems). 

Inspired by the PBEsol work, three of the present 
authors explored some other possible choices of con- 
straints for obtaining (3 and £i, 25 In one of these, 
PBE(G C , G x , LO), (3 and /i are both determined from gra- 
dient expansions, thus guaranteeing that this expansion 
is recovered, to the extent possible within the functional 
form of PBE, for both exchange and correlation. In an- 
other, PBE( J s , J r , LO), (3 and [i are both determined 
from jellium: /i from the jellium response function, as 
in PBE, and j3 from the jellium surface energy, as in 
PBEsol. Finally, PBE( J r , G x , LO) takes f3 from the jel- 
lium response function and /x from the gradient expan- 
sion of E x . The corresponding values of the parameters 
in each member of the PBE(/3, fi, LO) family are recorded 
in Table HI Additional information is given in Table I of 
Ref. [H. 



A priori one might expect that functionals such as 
PBE(G C , G Xl LO) and PBE( J s , J r , LO) that take (3 and fi 
from the same type of source, have the potential to ben- 
efit from error cancellation between the exchange and 
the correlation functional to a larger extent than func- 
tionals such as PBE and PBEsol that take them from 
different types of source. Also, one might anticipate that 
PBE(J S , J r , LO) should be rather good for simple met- 
als, as its takes both of its parameters from jellium, the 
paradigmatic model of such metals. These expectations 
were put to the test in Ref. [H, for atoms, molecules 
and solids. For each class of systems, a different rank- 
ing of functionals was found. Here we only record that 
for solids, where the calculations were done with the 
Siesta code^ original PBE performed worst of all. As 
expected, PBEsol provided significant improvement on 
PBE, but in spite of its name and the rationale behind 
its construction, it did not consistently provide the best 
performance for solids. Rather, best lattice constants 
were obtained from PBE(G C , G x , LO). It was not clear, 
however, to which extent this conclusion was affected by 
the pseudopotential approximation and the special basis 
functions employed in the Siesta code. 

In a second, independent, line of thought, the role of 
the Lieb-Oxford bound in functional construction has re- 
cently been reconsidered. Initial numerical and analyt- 
ical evidence^— suggested that the Lieb-Oxford esti- 
mate Alo = 2.273 could be tightened to a value close to 
A 2. Later, general arguments were given^ - that for 
three-dimensional systems this value should actually be 
Ael = 1-9555, where the subscript EL indicates that this 
is the exact value in the low-density limit of the electron 
liquid. This reduced value of A implies a corresponding 
reduction of n to 0.552. In our present notation, the 
resulting functional is denoted PBE(G C , J r , EL), and dif- 
fers from original PBE only in the value of A (or, equiva- 
lently, k) . This functional has been tested for a variety of 
molecular systems 2 ^ and it was found that PBE is rather 
insensitive to changes in A for covalently and ionically 
bound small molecules, a reduced, and thus, in principle, 
better, value of A producing slightly worsened energies 
and slightly improved bond lengths. 

In the present work we now tie up various open 
ends from these previous investigations, by implement- 
ing all ten functionals that can be obtained from the 
above-described combinations of /?, fj, and A, i.e., the 
complete family PBE(/3, /i, A), in the all-electron code 
WIEN2k^ and testing them systematically for a large 
set of 60 solids, comprising metals, semiconductors 
and insulators. 15 The PBE(/3,/i, A) functionals were also 
tested on a set of six molecules for the atomization en- 
ergy. This test set (called AE6) was proposed by Lynch 
and Truhlar— as a representative set of a much larger 
set of molecules. The molecules in the AE6 set are SiH-i, 
SiO, S 2 , C 3 H 4 , C 2 H 2 2 , and C 4 H 8 . 

The calculations on solids were performed with the 
WIEN2K code 2 - which solves the KS equations using the 
full-potential (linearized) augmented plane-wave and lo- 
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FIG. 1: (Color online) Relative error in the lattice constants of 28 elemental solids, obtained from LDA, original PBE = 
PBE(G C , J r , LO) and three versions of PBE that differ from it in just one constraint each, as described in the main text. Inset: 
Mean absolute relative error (mare) of the five functionals in the figure on this set of elemental solids. The Strukturbericht 
symbols (in parenthesis) are used for the structure: Al = fee, A2 = bec, and A4 = diamond. 



cal orbitals [FP-(L) APW+lo] method^, Because the FP- 
(L)APW+lo method is one of the most accurate methods 
to solve the KS equations it represents a good choice for 
testing xc functionals. The error in a calculated ground- 
state property is solely due to the approximate func- 
tional if good convergence parameters have been used. 
All calculations have been converged with respect to the 
number of k-points and the size of the basis set. Spin- 
orbit coupling for solids containing Ba, Ce, Hf, Ta, W, Ir, 
Pt, Au, Pb, and Th atoms has been taken into account. 
The experimental lattice constants are taken from Ref.fTol 
and are corrected for zero-point anharmonic expansion. 
The calculations on molecules were done with the de- 
Mon code22 which uses Gaussian basis sets. The very 
large uncontracted basis sets developed by Partridg o 40 ' 41 
were used. 

The statistical quantities that will be used for the anal- 
ysis are displayed below, where p£ alc and p^ xp are the cal- 
culated and experimental values of the considered prop- 
erty (either the lattice constant or the atomization en- 
ergy) of the ith solid or molecule of the testing set: 

mean error (in A or in kcal/mol), 



■E 



„calc 



(3) 



the mean absolute error (in A or in kcal/mol), 
1 " 

mae = - £ k alc ~ Pi 



cxp I 



(4) 



cxp 



the mean relative error (in %) 

1 n calc 

mre = - V 100^ 
and the mean absolute relative error (in %), 



1 

mare = — >^ 100 
n ^ 

i=l 



~calc 



cxp 



cxp 



(5) 



(6) 



The spread (in %), defined as 



spread = max 100 



calc 



cxp 



- min 100 



calc 



cxp 



Pi 



(7) 



will also be discussed. The smaller the spread, the more 
predictable a functional behaves. A large spread, by con- 
trast, indicates a more erratic behaviour. In situations 
where a single bad value can be problematic, it may be 
wiser to choose a functional with a small spread than one 
with a low mean error, as the latter may still be way off 
in isolated cases. 



III. ANALYSIS OF SINGLE-CONSTRAINT 
CHANGES WITH RESPECT TO PBE 

The calculated lattice constants for all functionals con- 
sidered in this work are given in Table SI of the supple- 
mentary EPAPS material4? Graphical representations of 
these results are also given in Figs. SI, S2, and S3. 
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In a first step, we focus our analysis on a subset of 
functionals that differ from original PBE in the choice 
of just one constraint. Specifically, these functionals are 
PBE( J s , J r , LO), differing in the constraint for determin- 
ing j3, PBE(G C , G x , LO), differing in the constraint cho- 
sen for fi, and PBE(G C , J r , EL), differing in the choice of 
A. Initial focus on just these functionals is useful because 
it allows to separately assess the influence of each change 
relative to PBE. 

Figures[T]and[2]show the signed relative errors for all 60 
solids, for the three functionals just described as well as 
for original PBE and LDA. While original PBE has the 
same (modest) performance for elemental solids as for 
compounds, all other functionals work better for com- 
pounds. The difference is particularly pronounced for 
LDA and PBE(G C , J r , EL), which perform significantly 
better for compounds than for elemental solids. In both 
classes, however, PBE(G C , G x , LO) achieves the lowest 
mare among this subset of functionals. 

Figures [T]and [5] also show very clearly the known trend 
of LDA to underestimate the lattice constants (negative 
relative errors) and of PBE to overestimate the lattice 
constants (positive relative errors). The most interest- 
ing fact, which is not at all obvious from the way the 
various functionals were constructed, is that all changes 
of parameters relative to PBE produce significantly better 
lattice constants. (This remains true even if all the other 
possible combinations are included.) Since the nature 
and source of each modified parameter are completely 
different in each of the three cases, this seems to indicate 
that the original PBE choice was rather unfortunate for 
lattice constants, as reasonable changes to any of its pa- 
rameters end up improving the results. 

Comparing the absolute size of the change resulting 
from each modified parameter, we immediately conclude 
from Figs. [T]and[2]that PBE is most sensitive to changes 
in \i and least sensitive to changes in \M The relative 
impact of f3 and \i is quite reasonable, because f3 appears 
in the correlation energy and fi in the exchange energy. 
Since the exchange energy in ordinary solids is larger than 
the correlation energy the result should indeed depend 
more sensitively on changes of that quantity. 

Regarding changes only in A, we see that a change from 
the Lieb- Oxford value to the electron-gas value has, for 
most solids, the smallest effect on PBE of the three tested 
parameter changes. This is consistent with what was pre- 
viously observed for molecules^, The exception is for the 
alkali metals, where this change has the biggest effect. 
The reasons for this behaviour becomes clear from the 
analysis given by Haas et al. (Ref. l24h . In this paper an 
"important region" was defined which is to a large extent 
responsible for the changes in lattice parameters of dif- 
ferent functionals. For closed packed solids (like elements 
in the fee or bee structure) this region is the separation 
between the outermost core ("semi-core") and the va- 
lence states and for the alkali metals the reduced density 
gradient s = |Vp| / ^2 (37r 2 ) p 4 / 3 ^ in this "important 
region" is much larger (even above s — 2 for Li) than for 



other elements (e.g., in bee V, s max = 0.9). Obviously, 
a change in A modifies the enhancement factor F xc (r s , s) 
much more for large s, while a change in /x influences 
F xc (r s , s) predominantly in the low s region. 



IV. FULL STATISTICS FOR ALL TEN 

PBE(^ lA t,A) FUNCTIONALS 

In this section we present a statistical analysis of all ten 
functionals of the PBE(/3, n, A) family, as well as of LDA. 
We do not include numerical data on other GGA-type 
functionals, such as SOGGA 22 and WCS and neither on 
alternative functionals such as AMOS^ and TPSS meta- 
GGA^ as these were already investigated on the same 
test set in Ref. [HI, and on a smaller set (using different 
codes and basis sets) in Refs. ll8ll35ll36L However, in our 
discussion of global trends we compare with results and 
conclusions from those references. We do not consider 
earlier variations of PBE, such as revPBE (proposed in 
the Comment of Zhang and Yang£) and RPBE 4 ^ which 
were not designed for solids and tend to worsen PBE for 
extended systems (see, e.g., Ref. |H for RPBE). 

The choice of the best performing functional for lat- 
tice constants depends on the measure of error selected 
to judge the performance of the functionals. In terms 
of the mean error, PBE( J r , G x , LO) achieves a spec- 
tacularly low error of -0.002 A, followed by PBEsol 
and PBE( J s , J r , EL). These same three functionals also 
achieve the best mean relative errors. However, such 
low errors can in part be due to the result of error 
cancellation in taking the averages. The mean abso- 
lute error, which is not influenced by error cancellation, 
favours PBEsol, closely followed by PBE(G C , G x , LO) 
and PBE( J r , G x , LO). The same three functionals also 
achieve the lowest mean absolute relative error. 

In view of the very small differences in the mares of 
some of the functionals, one must ask how significant 
these differences are, considering both, the numerical ac- 
curacy of the theoretical results and the accuracy of the 
experimental data, where not for all cases high quality 
low-temperature lattice parameters and good zero-point 
energy corrections are available. In fact, changes in mare 
of 0.05% are about the limit of theoretical accuracy, i.e., 
we expect to have an absolute precision of about 0.005 
A. Therefore, we consider PBEsol, PBE( J r , G x , LO) and 
PBE(G C , G x , LO) to perform equally well in terms of the 
mare. 

Three of the four functionals with good mre or mare 
take the original Alo and A* from the gradient expansion 
for exchange (G x ), but differ in the parameter of the 
correlation functional, f3. As the numbers show, the value 
of f3 turns out to be almost irrelevant, while the value of 
/j appears to be responsible for the improved behaviour. 
The gradient expansion for exchange is thus seen to be 
the key ingredient in functionals that deliver good lattice 
constants. When /i is taken a bit larger (e.g., J r ), both 
(3 and A must be set to small values (f3 — J s and A = 
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FIG. 2: (Color online) Relative error in the lattice constants of 32 compound solids, obtained from LDA, original PBE = 
PBE(G C , J r , LO) and three versions of PBE that differ from it in just one constraint each, as described in the main text. Inset: 
Mean absolute relative error (mare) of the five functionals in the figure on this set of compound solids. The Strukturbericht 
symbols (in parenthesis) are used for the structure: Bl = rock-salt, B2 = cesium-chloride, B3 = zinc-blende, and CI = fluorite. 



Ael) to obtain a similar performance [PBE( J s , J r , EL)]. 
This observation, in retrospect, vindicates the PBEsol 
approach^! and attests to the solidity of the basic insight 
presented there regarding the relevance of the gradient 
expansion of the exchange energy for solids. 

In terms of the spread [Eq. (J7J)], the best performer is 
PBE(J S , J r ,LO), followed by PBE(J r , G x , LO). Unfor- 
tunately, the best performers with regard to mean errors 
and with regard to the spread are not always the same. 
In particular PBE(G C , G x , LO) has a rather large spread, 
but PBE( J r , G x , LO) appears to be a reasonable compro- 
mise, doing well according to all three criteria. 

Next, we compare our present conclusions to several 
different sets of earlier calculations, testing specific mem- 
bers of the full family, employing other test sets, or other 
implementations and basis functions. 

In Ref. HE], three of us tested the five functionals 
PBE(/3, \i, LO) on a set of 13 solids. These functionals 
were implemented in the Siesta code^ This code, by 
design, aims at the electronic structure of very large sys- 
tems, where all-electron calculations, even with simple 
functionals, would be prohibitively expensive. To this 
end, it makes use of specially designed localized numer- 
ical basis functions, and pseudopotentials. As a conse- 
quence, it does not attain the same high accuracy as all- 
electron codes, such as WIEN2k^ and the absolute size 
of the errors is larger for Siesta than for WIEN2k. 

Nevertheless, the resulting error statistics is rather 
similar (although not identical). In particular, both 
the Siesta and the WIEN2k calculations identify orig- 
inal PBE and PBE( J s , J r , LO) as the worst performers 
for lattice constants of all PBE(/3, /i, LO) functionals, and 



PBEsol, PBE( J r , Gz, LO), and PBE(G C , G x , LO) tied as 
the best. Among these best performing functionals, 
Siesta and WIEN2k produce a different ranking, with 
Siesta preferring PBE(G C , G x , LO), which according to 
WIEN2k is beaten by a small margin by PBEsol and 
PBE{J r ,G x ,LO). 

Interestingly, all-electron calculations for solids per- 
formed with the Gaussian code^S also indicate that 
PBE(G C , G x , LO) produces better lattice constants than 
PBEsol (see Table SIV in the supplementary EPAPS ma- 
terial of Ref. [U) which is in line with the Siesta re- 
sults. However since neither pseudopotentials (Siesta) 
nor Gaussian basis functions (Gaussian) are as accurate 
for solids as all-electron FP-(L)APW+lo calculations, 
and since all differences are rather small, we still regard 
these functionals as essentially tied. 

In Ref. H^, two of us with Samuel B. Trickey tested 
the functional PBE(G C , J r ,EL), which differs from orig- 
inal PBE only in the reduction of A, corresponding to a 
tighter (and thus presumably better) Lieb-Oxford bound. 
The calculations were done for a set of small molecules. 
The reduction of A was found to slightly improve inter- 
atomic distances. In parallel, from Siesta pseudopoten- 
tial calculations^ for the 13 solids of Ref. [25| we found 
that the same improvement occurs also for the other 
members of the PBE(/3, /i, A) family, all of which pro- 
duce better lattice constants when a tighter Lieb-Oxford 
bound is enforced. By contrast, in the WIEN2k calcu- 
lations only the badly performing functionals [original 
PBE and PBE( J s , J r , LO)] benefit from a reduced value 
of A, while the mare of the other three functionals grows 
if A is reduced. Consistently with what was speculated 
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TABLE I: Parameters of the eleven functionals under investigation in this work, and a statistical summary [Eqs. l|3" |l -l|? |) ] of 
their performance for lattice constants of 60 solids. Lattice constants differing by less than ~ 0.005 A and mares differing by 
less than « 0.05% should be considered equivalent. The lattice constants are given in Table SI of the supplementary EPAPS 
material.— 



Functional 



P 



/' 



A' 



(A) mae (A) 



(%) 



(%) spread (%) 



LDA 

PBE(G C , J r ,LO) a 
PBE(J s ,G I ,LO) 6 
PBE(J 3 , J r ,LO) c 
PBE(G c ,G x ,LO) c 
PBE(J r ,G ;c ,LO) <; 
PBE(G C , Jr,EL) d 
PBE(J S ,G„EL) 
PBE(J 3 , J r , EL) 
PBE(G C , Gz, EL) 
PBE(J r ,G ;c ,EL) 



-0.060 0.060 -1.37 1.37 4.95 

0.067 0.21951 2.273 0.049 0.053 1.02 1.14 3.68 

0.046 0.12346 2.273 -0.007 0.028 -0.21 0.64 3.71 

0.046 0.15133 2.273 0.014 0.033 0.25 0.72 3.41 

0.067 0.12346 2.273 -0.012 0.029 -0.31 0.67 4.30 

0.038 0.12346 2.273 -0.002 0.030 -0.09 0.67 3.56 

0.067 0.21951 1.9555 0.014 0.037 0.28 0.80 3.73 

0.046 0.12346 1.9555 -0.023 0.033 -0.55 0.76 4.00 

0.046 0.15133 1.9555 -0.008 0.032 -0.20 0.71 3.87 

0.067 0.12346 1.9555 -0.028 0.034 -0.67 0.78 3.75 

0.038 0.12346 1.9555 -0.018 0.034 -0.44 0.75 4.08 



a This is PBE.1 
<This is PBEsol^l 

c These are the three functionals proposed in Ref. 25. 
d This is the functional of Ref. [H. 

e A = 2.273 and 1.9555 correspond to K = 0.804 and 0.552, respec- 
tively [see Eq. J2}]. 



in Ref. [29|, this indicates that the algebraic form of PBE 
is too restricted to systematically benefit from a tighter 
bound. 

In Ref. [H], three of us employed the same set of 
60 solids to assess the performance of the AM05, WC, 
PBEsol, SOGGA, and the meta-GGA TPSS functionals 
compared to the older LDA and PBE. From this com- 
parison PBEsol emerged as the functional with the low- 
est mare over all 60 solids, tied with the WC, closely 
followed by SOGGA and AM05, and more distantly by 
TPSS, PBE and LDA, in this order. WC, SOGGA, and 
AM05 turn out to have mares in the same range as the 
members of the PBE(/3, /i, A) family (including PBEsol) 
although they differ from the original PBE by more than 
just the values of parameters. 

The hybrid functional B3LYP ) 10 ' 11 which is very pop- 
ular in quantum chemistry, was showr*^ to overestimate 
lattice constants by at least as much as PBE, and is thus 
not competitive with any of the functionals under study 
here. 

Interestingly, SOGGA turns out to be a very good 
functional making use of the tighter Lieb-Oxford bound, 
using Ael instead of Alo- SOGGA achieves a lower mare 
(0.68%) than any of the five functionals PBE(/3, /i, EL), 
whose mare ranges from 0.71 to 0.80%, but unfortunately 
its mre is twice as large as that of PBE( J s , J r , EL). This 
indicates, one more time, that the functional form of PBE 
must be changed to fully benefit from a tighter Lieb- 
Oxford bound, and hints that the form of SOGGA may 
be a suitable starting point for this purpose. 



V. ANALYSIS FOR CLASSES OF SYSTEMS 

The above discussion was based on the statistical data 
of the full set of 60 solids. However, it can also be in- 
teresting to analyze the results for a particular class of 
solids, therefore, below we discuss the performance of all 
11 functionals separately for certain classes of solids, for 
which Table SI and Fig. SI of the supplementary EPAPS 
material^ are useful. 



A. Elemental solids 

Let us start out the discussion with the alkali met- 
als. For a given /3 and /i, the reduction of A from 
LO to EL leads to significantly smaller lattice constants 
aQ. As mentioned before, this effect is particularly 
strong for the alkali metals and thus all PBE(/3, /i, A) 
functionals with a tighter Lieb-Oxford bound underes- 
timate the lattice constant. Actually, from Fig. [1] and 
SI we can see that for the alkali metals (and also the 
alkali-earth metals and the compounds with these el- 
ements) the difference between the PBE (A = Alo) 
and LDA (A = 1) relative errors is large, an effect 
which is (at least partially) due to the large values of 
s (which make A important) in the region of separation 
between semi-core and valence electrons^ By compar- 
ing the results with fixed /i = G x and a variation of 
/? [PBE(J s ,G^,LO/EL) with PBE(G C , G x , LO/EL) and 
PBE( J r , G x , LO/EL)] we note that an increase of (3 from 
J s to G c increases the lattice constant more than a re- 
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duction of J r decreases it. J r worsens the underestima- 
tion of ao from PBE( J s , G x , LO/EL) while G c reduces 
it. This trend, namely that a (PBE(G c , G x , LO/EL)) 
> <zo(PBE(J s , Gjb, LO/EL)) > a (PBE(J r , G Xl LO/EL)) 
can also be observed for most group IIA elements, but in 
all other cases the trend is inverted. On the other hand, 
when fi is taken as J r [and therefore depends on /3, see 
Eq. ((IJ)], an increase of — J s to (3 = G c [comparing 
PBE(J S , J r , LO/EL) with PBE(G C , J r , LO/EL)] leads to 
increased lattice parameters for all classes of compounds, 
not just the group IA and IIA elements. This effect is 
most pronounced in K and Rb and reduces the large over- 
estimation of PBE(G C , J r , LO) for these two compounds. 
Usually, a reduction of fj, = J r to /j, = G x (at any (3 
and A) [PBE(G C , J r , LO/EL) and PBE(G C , G x , LO/EL)] 
leads to much smaller lattice constants, but for group IA 
elements this general trend is not true for Li (and in very 
few cases for Na and K) . In general, for group IA elements 
PBE( J s , G x , LO) is the most accurate PBE(/3, ^, A) func- 
tional followed by PBE(G C , G x , LO). 

For the elements of group IIA the original PBE 
[PBE(G C , J r , LO)] gives already quite satisfactory results 
and thus all modifications of the PBE(/3, /x, A) functionals 
lead to strong underestimations of the lattice constants. 
Tightening the Lieb-Oxford bound increases the absolute 
relative error by about 1.5%. Reduction of /x has an even 
larger negative effect. Changing (i (at hxed [i = G x ) 
has a fairly small effect (in particular for Ba) and a non- 
uniform trend. 

The lattice parameters of group IVA elements are very 
well described by standard LDA (only for Pb there is 
a significant underestimation), while the original PBE 
functional yields ao almost 3% too large for Sn and Pb. 
In addition, PBE shows a strong tendency for larger over- 
estimation of ao for heavier elements. Since LDA is so 
good, only the "weakest" GGAs, i.e., where the exchange 
parameter fi is small (G x ), can compete with LDA. 
When in addition also the reduced Ael is used and/or 
(3 = G c is kept large, functionals like PBE(G C , G x , EL), 
PBE(J S ,G X ,EL) or PBE(G C , G x , LO) (in that order) 
perform very well. 

The trend for the 3d transition metals (TM) is very 
similar to that for the group IIA elements. Since already 
original PBE is rather accurate for the 3d TM (except for 
Cu), no overall improvement can be expected when one 
(or several) of the parameters is reduced. A tighter Lieb- 
Oxford bound [PBE(G C , J r , EL)] improves the situation 
for Cu, but worsens all other cases. A reduction of // has 
an even larger negative effect, but PBE( J s , J r , LO) is the 
best modified A = Alo functional. None of the modifi- 
cations seems to be able to break the trend that lattice 
parameters of early 3d TM are even more underestimated 
than of late ones. 

The lattice constants of the 4d transition metals are 
overestimated by original PBE (the overestimation in- 
creases for later TM) and slightly underestimated by 
PBEsol (getting more accurate with increasing nuclear 
charge). Using a tighter Lieb-Oxford bound (Ael) there- 



fore improves PBE (worsens PBEsol), but this modifica- 
tion alone is not enough. An additional reduction of /3 
from G c to J s [PBE( J s , J r , EL), which also reduces the ef- 
fective /j] leads to pretty accurate results, while the (3 re- 
duction alone [PBE(J S , J r ,LO)] is not sufficient. Similar 
good results can also be reached when both, \x and /3 are 
reduced to G x and J r , respectively [PBE( J r , G x , LO)]. 
Most interestingly, these two modifications can also sig- 
nificantly reduce the trend towards larger lattice param- 
eters for later TM and are thus an improvement for all 
4d elements. Any further combination with reduced /3, /z, 
and A underestimates the lattice parameters. 

For the 5d TM, original PBE overestimates ao and 
for the latest 5d element (Au), the error reaches more 
than 2%. The best functional for the 5d elements is 
PBE(G C , G x , LO), where only the exchange factor [i is 
strongly reduced, but (i (and A) are kept at the large 
values. As for the 4d series, the trend towards larger 
lattice parameters for late TM elements is more or less 
completely broken. Similar, but slightly overestimated 
ao can be obtained when both, (3 and A are also reduced 
[PBE( J s , G x , EL)]. Reduction of A alone or intermediate 
values for fi still overestimate the lattice parameters. 

For the heaviest element of our testing set, the 5/ el- 
ement Th, the original PBE gives the best result (still 
underestimating ao slightly), while, e.g., PBEsol leads to 
a more than 2% too small lattice parameter. 



B. Compounds 

Most prior discussed trends (Sec. IV A[) can to some 
extent also be observed for compounds. Sometimes the 
combined effect of two elements may lead to some kind of 
cancellation, or in other cases, one element may dominate 
the effect. We will discuss below the effects starting with 
the very ionic group I- VII and II- VI compounds, then 
covering the more covalently bound group III-V and TM- 
compounds. 

For the ionic compounds it is obvious that the anion 
(the tails of the valence p electron density) plays the 
major role in determining the lattice parameter. This 
was shown in Ref. HU, but is also obvious by compar- 
ing the lattice parameters for e.g. metallic Li and LiF. 
Using PBE there is a small underestimation of ao for 
Li, but the lattice parameter of LiF is too large by al- 
most 3%. In addition, the anion changes dramatically 
the results: For fluorides a large overestimation is ob- 
tained for all PBE(/3, fi, A) functionals, while for chlo- 
rides (and even more for bromides, not included here) 
this behaviour is corrected or some underestimation can 
be found. The change to a tighter Lieb-Oxford bound 
has a rather strong effect (often stronger than a re- 
duction of /i) and reduces most errors of the IA-VIIA 
compounds. While nearly all PBE(/3, fi, LO) functionals 
overestimate the lattice constants, a tighter Lieb-Oxford 
bound may lead to small underestimations for some chlo- 
rides. Nevertheless PBE( J r , G x , EL) is the best perform- 
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TABLE II: Parameters of the eleven functionals under investigation in this work, and a statistical summary [Eqs. ©-(O] of 
their performance for the atomization energy of the set AE6 of six molecules.^ The atomization energies are given in Table 
SII of the supplementary EPAPS material.— 



Functional /3 fj, \ e me (kcal/mol) mae (kcal/mol) mre (%) mare (%) 



LDA 


76.3 


76.3 


16.9 


16.9 


PBE(G C , J, ,LO) a 0.067 0.21951 2.273 


12.0 


15.1 


3.4 


4.4 


PBE{J 3 ,G x ,L0) b 0.046 0.12346 2.273 


35.1 


35.1 


8.3 


8.3 


PBE(J 3 , J r ,LO) c 0.046 0.15133 2.273 


28.5 


28.7 


6.9 


6.9 


PBE(G c ,G x ,LO) c 0.067 0.12346 2.273 


31.0 


32.7 


7.6 


8.2 


PBE(J r ,G x ,LO) c 0.038 0.12346 2.273 


36.4 


36.4 


8.5 


8.5 


PBE(G C , Jr,EL) d 0.067 0.21951 1.9555 


26.6 


28.5 


6.5 


7.1 


PBE(J S ,G X ,EL) 0.046 0.12346 1.9555 


43.5 


43.5 


10.1 


10.1 


PBE(J 3 , J,-, EL) 0.046 0.15133 1.9555 


38.9 


38.9 


9.1 


9.1 


PBE(G C , Gx, EL) 0.067 0.12346 1.9555 


39.4 


40.4 


9.4 


9.7 


PBE(J r , Gz, EL) 0.038 0.12346 1.9555 


44.8 


44.8 


10.3 


10.3 



"This is PBEX 
6 This is PBEsol^i 

These are the three functionals proposed in Ref. 25. 
d This is the functional of Ref. 29. 

e A = 2.273 and 1.9555 correspond to K = 0.804 and 0.552, respec- 
tively [see Eq. Idl- 



ing PBE(/3, zx, A) functional for this class of compounds. 

For the IIA-VIA compounds standard PBE overes- 
timates the lattice constants. A reduction from f3 = 
G c to (3 — J s at xx = J r [comparing PBE with 
PBE(J S , J r , LO/EL)] improves PBE significantly, be- 
cause the change of (3 reduces the effective xx and there- 
fore PBE( J s , J r , LO) becomes the most accurate func- 
tional for this group of compounds. The reduction from 
/x = J r to /x = G x at (3 — G c [e.g., from PBE to 
PBE(G C , G x , LO)] also lowers the lattice constants signif- 
icantly leading to quite well performing functionals with 
Alo , but with Ael the correction overshoots and leads 
to some underestimation for MgS and CaO. For a fixed /x 
the variation of f3 hardly modifies the results. Only MgO 
[which is overestimated by all PBE(/3,/x, A) functionals] 
benefits from a tighter Lieb-Oxford bound in all cases. 

All semiconducting IIIA-VA compounds have signifi- 
cantly too large lattice parameters with PBE. A tighter 
Lieb-Oxford bound is advantageous, but the effect alone 
is too small and a strong reduction of /x (to G x ) is es- 
sential. For /x = G x the parameter for correlation has a 
fairly large effect (in particular for the heavier elements) 
and with (3 = G c the PBE(G C , G x , LO) functional is very 
accurate for all semiconductors. 

The metallic transition metal compounds (mainly car- 
bides and nitrides, but also three intermetallic com- 
pounds) are fairly well described by standard PBE and 
the slight overestimation of lattice constants can be 
reduced by weak modifications. The best functionals 
are obtained either by switching A to the tighter EL 
limit [PBE(G C , J r , EL)], or by a modest reduction of (3 
[PBE(J S , J r , LO)] (probably because this reduces also the 
effective /x for /i = J s ). A stronger reduction of /x or a 



combination of /x and A reductions leads to too small 
lattice parameters. 



VI. ATOMIZATION ENERGY OF MOLECULES 

In this section we present the performance of the 
PBE(/3, /x, A) functionals on the atomization energies of 
molecules using the representative AE6 test seti^ Table 
SII of the supplementary EPAPS material^ gives the 
calculated values and Table |TT] gives a summary of the 
corresponding statistical errors. The me and mae (mre 
and mare) quantities are very similar since all function- 
als (maybe except original PBE) always overestimate the 
atomization energy. We can see that PBE(G C , J r , LO) 
(original PBE) is the best and PBE( J ri G Xl EL) is the 
worst functional of the PBE(/3, /x, A) family for the at- 
omization energy of molecules. Switching A to a tighter 
bound (A = Ael) has a rather strong degrading effect, 
probably because the atomization energies depend a lot 
on regions in space with large effective density gradi- 
ent ,s. In Fig. |3] we compare the mre and mare for 
solids (lattice constants) and molecules (atomization en- 
ergies) versus the PBE(/3,/i, A) functionals. As expected 
the mre behavior of solids is opposite to that of the 
molecules. Functionals leading to larger lattice constants 
(larger overestimation) lead to smaller atomization en- 
ergies (smaller overestimation) and no functional of the 
current PBE(/3, ix, A) family can describe both quantities 
in a satisfying way. 
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FIG. 3: Mean relative error and mean absolute relative error 
in the lattice constants (upper panel) and atomization ener- 
gies (lower panel) of the PBE(/3, fi, A) functionals. 



VII. CONCLUSIONS 

From all of the above we conclude that to obtain pre- 
cise lattice constants of solids it is not necessary (and 
in some cases even detrimental) to switch from the PBE 
family of functionals (differing from original PBE only 
through the choice of parameters) to functionals that also 
differ in the form of the exchange enhancement factor 
(SOGGA and WC) or that employ different design prin- 
ciples (AM05) or further ingredients (meta-GGA TPSS). 



These more complex functionals have many merits and 
interesting features, but apparently these features are not 
required to produce very accurate lattice constants. 

In fact, even the simplest possible modification of orig- 
inal PBE, the change of one single parameter (taking /i 
from the gradient expansion for exchange instead of from 
the jellium response function) already produces a func- 
tional whose lattice constants are, to within the error 
bars of the WIEN2k code, as good or better than those 
of any of the other tested functionals: PBE(G C , G x , LO). 
(As pointed out above, Siesta and Gaussian calculations 
sustain this claim.) 

A change of two parameters relative to PBE produces 
PBE(J r ,G x ,LO) and PBEsol = PBE( J s , G x , LO), the 
former having the same mare but a lower spread, rela- 
tive to PBE(G C , G x , LO); and the latter a still slightly 
lower mare (although the improvement is smaller than 
our estimated error bar) at the expense of a slightly larger 
spread. Any of these three functionals can be recom- 
mended as a useful and reliable GGA for lattice con- 
stants of solids, requiring only minimal changes to exist- 
ing implementations of PBE and attaining much higher 
accuracy than PBE and LDA, and also than many more 
complex functionals. 

On the other hand, any modification of the original 
PBE functional which improves lattice parameters of 
solids, increases significantly the error in the atomization 
energy of molecules. The overestimation of this quan- 
tity by PBE of mae = 15.1 kcal/mol becomes 2—3 times 
larger with the modified functionals and we must con- 
clude that (at least) GGAs of PBE-form cannot describe 
well lattice parameters of solids and atomization energies 
of molecules simultaneously (see also Refs. [22j andl48l). 
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